Proyectando el comportamiento de la soja¶
Competencia Metadata 2019
https://metadata.fundacionsadosky.org.ar/competition/11/
Objetivos¶
- Lograr el mejor fit para la serie de retornos simples diarios $\frac{p_1 - p_0}{p_0}$. Serie desde la Fecha de Cierre del Concurso (FCC 29/9/2019) + 10 días hábiles.
- Lograr la mejor proyección para el último valor de la serie contínua de la soja al fin del concurso. Cierra a la FCC + 10 días hábiles.
- Lograr el mejor fit para la serie de retornos simples diarios de 4 semanas. Desde FCC - 10 días hábiles, hasta FCC + 10 días hábiles.
Dataset oficial¶
Dataset la información que contiene este archivo es hasta el día 29/8/2019. Por ende para calcular el score circunstancial se entiende que usted está proyectando hasta el día 12/9/2019. Esta información se irá actualizando todas las semanas hasta la última actualización que será el día 27/9/2019.
Descripción¶
FechaFecha de referencia para el precio.OpenPrecio de apertura del día.HighPrecio máximo del día.LowPrecio mínimo del día.LastPrecio último operado del día.CierrePrecio de ajuste del día. Ésta es la serie a proyectarse.Aj.Dif.Diferencia nominal respecto del día anterior.MonMoneda de denominación del contrato.Oi.VolInterés abierto del contrato.Oi.Dif.Diferencia del interés abierto respecto del día anterior.Vol.Ope.Volumen Operado medido en contratos.UnidadUnidad en que se miden los contratos.DolarB.N.Precio del dólar del Banco de la Nación Argentina.DolarItauPrecio del dólar del Banco Itaú.Diff.SemDiferencia Semanal.
Formato de respuesta¶
Debe enviarse un archivo en formato csv sin encabezado con 4 columnas y 20 filas.
- Primer columna debe contener al número de fila (int).
- Segunda columna debe contener las fechas correspondientes en formato "dd/mm/YYYY."
- Tercer columna debe contener los retornos (float).
- Cuarta columna debe contener el precio (cierre del día) del contrato (float).
Las filas 1 a 10 corresponden a FCC - 10 días hábiles y las filas 11 a 20 corresponden a la proyección futura.
Exploración¶
%config InlineBackend.figure_format = 'retina'
import os
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from pandas.plotting import register_matplotlib_converters
import pyflux as pf
import seaborn as sns
import statsmodels.api as sm
warnings.simplefilter('ignore')
register_matplotlib_converters()
sns.set(rc={'figure.figsize': (12,8)})
data_dir = 'data'
date_parser = lambda x: pd.datetime.strptime(x, "%d/%m/%Y 12:00:00 a.m.")
df = pd.read_csv(os.path.join(data_dir, 'datasetRofex2.csv'),
parse_dates=['Fecha'],
index_col='Fecha',
date_parser=date_parser)
df.head()
Empezamos con un plot del precio de cierre en cada día.
df['Cierre'].plot(title='Precio diario de cierre de Soja FAS');
El gráfico sugiere que se trata de una serie no estacionaria, con un trend marcado. Podemos confirmar esto con un gáfico de autocorrelación. Más aún, con el plot de autocorrelacion parcial, podemos ver que la maxima autocorrelacion esta en el primer lag.
fig, (ax1, ax2) = plt.subplots(2)
sm.graphics.tsa.plot_acf(df['Cierre'], ax=ax1)
sm.graphics.tsa.plot_pacf(df['Cierre'], ax=ax2);
Si la serie fuera estacionaria, uno esperaría que la autocorrelación decreciera rápidamente al aumentar el lag, pero esto no es así. Nuestro objetivo, sin embargo, no es predecir el precio de cierre sino el retorno diario, que es (casi) la serie diferenciada. Veamos la autocorrelación de estos retornos.
df['retorno'] = df['Cierre'].pct_change()
retornos = df['retorno'].dropna()
df['retorno'].plot(kind='hist', bins=50, title='Frecuencia de retornos simples diarios', density=True);
df['retorno'].dropna().plot(title='Retornos diarios');
Repetimos los plots de autocorrelación para la serie de retornos diarios.
fig, (ax1, ax2) = plt.subplots(2)
sm.graphics.tsa.plot_acf(df['retorno'].dropna(), ax=ax1)
sm.graphics.tsa.plot_pacf(df['retorno'].dropna(), ax=ax2);
Se ve claramente que al tomar retornos, la autocorrelación baja muy rápido al aumentar el lag, lo que nos dice que estamos en presencia de un proceso más estacionario.
Preliminares¶
Para determinar el grado de estacionaridad de la serie, calculamos el exponente de Hurst, una medida que caracteriza el grado de memoria de la serie.
La idea es considerar la varianza de la serie de lags. Para un lag de orden $\tau$, tenemos:
El exponente de Hurst $H$ queda definido como:
$${\rm Var}(\tau) = \tau^{2H}$$De acuerdo al valor de dicho exponente, la serie en cuestión es:
- $H < 0.5$ - la serie revierte a la media
- $H = 0.5$ - la serie es un movimiento browniano geométrico
- $H > 0.5$ - la serie tiene una tendencia
Para una explicación mas detallada, consultar el siguiente post.
# Código adaptado de https://www.quantstart.com/articles/Basics-of-Statistical-Mean-Reversion-Testing
def hurst(ts):
lags = range(2, 100)
tau = [np.sqrt(np.std(np.subtract(ts[lag:], ts[:-lag]))) for lag in lags]
poly = np.polyfit(np.log(lags), np.log(tau), 1)
return poly[0] * 2.0
hurst(df['Cierre'])
hurst(df['retorno'])
Vemos que la serie de precios diarios se aproxima a un movimiento browniano, mientras que la serie de retornos presenta una reversión a la media.
Método ARIMA¶
Nuestra herramienta principal para la predicción de los retornos será el método ARIMA, que explicamos brevemente a continuación.
ARIMA (Autoregressive integrated moving average) es un modelo aplicado a series de tiempo que considera al valor actual de la serie como una regresión sobre valores previos y al error de regresión como una combinación lineal de errores previos. Más precisamente, si la serie la notamos $\{Y_t\}_{t \in \mathbb{N}}$, el valor a tiempo $t$ se modela por
$$Y_t = \phi_0 + \sum_{i=1}^{p}{\phi_i Y_{t-i}} + \sum_{i=1}^{q}{\theta_i \varepsilon_{t-i}} + \varepsilon_t$$donde los $\varepsilon_i$ con $i\neq t$ son errores de regresión, $\varepsilon_t$ es ruido blanco y los coeficientes $\phi_i$ y $\theta_i$ se ajustan típicamente a través de métodos estadísticos. La cantidad de pasos en el pasado que se consideran forma un par de hiperparámetros del modelo (los números $p$ y $q$), con un tercero (generalmente notado $d$) que determina cuántas veces diferenciar la serie antes de aplicarle lo anterior. Todo esto da un método $ARIMA(p,d,q)$. En nuestro caso tomamos $d=0$ porque aplicamos el método a los retornos, que ya son lo suficientemente estacionarios. Si lo quisiéramos aplicar al precio de cierre tomaríamos $d=1$ para diferenciarlo, obteniendo una serie muy similar a los retornos pero a otra escala (porque en cada paso no estaríamos diviendo por el precio del día anterior).
Una vez obtenido el fit por ARIMA, usamos el método GARCH para hacernos una idea de la volatilidad de los residuos y corroborar que ésta es (aproximadamente) constante.
GARCH (Generalized autoregressive conditional heteroskedasticity) es otro modelo aplicado a series de tiempo para describir la varianza del error entre el valor observado y el que se había previsto (en nuestro caso, el término $\varepsilon_t$ que se consideraba ruido blanco). Típicamente esta varianza se escribe como combinación lineal de sus valores previos, como en un método ARMA (la $i$ de ARIMA refiere al proceso de diferenciar la serie para volverla estacionaria).
Para decidir qué hiperparámetros de ARIMA tomar usamos el método auto_arima de la librería pmdarima, que ajusta el método para diferentes valores de $p$ y $q$, los compara usando el criterio de información de Akaike y el criterio de información bayesiano y arroja el mejor fit.
import pmdarima as pm
stepwise_fit = pm.auto_arima(retornos, start_p=1, start_q=1,
max_p=3, max_q=3,
start_P=0, seasonal=False,
d=None, trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True)
stepwise_fit.summary()
#ARIMA
p = 1
q = 1
#GARCH
r = 1
s = 1
#Número de predicciones a futuro
n = 13
model_arima = sm.tsa.statespace.SARIMAX(retornos.values,order=(p,0,q), seasonal=False)
fit = model_arima.fit()
mean_prediction = fit.forecast(n)
Graficamos las predicciones de ARIMA, primero en los retornos y después en el precio de cierre.
pd.Series(mean_prediction).plot();
pd.Series((1 + mean_prediction).cumprod() * df['Cierre'][-1]).plot();
residuos = fit.resid
model_garch = pf.GARCH(residuos,r,s)
x = model_garch.fit()
x.summary()
model_garch.plot_predict(h=10)
model_garch.plot_predict_is(h=50,figsize=(15,5))
Comparamos ahora las predicciones obtenidas con los valores reales del período 30/08 - 17/09, los 13 días hábiles inmediatamente posteriores a nuestra serie.
futuros = pd.read_csv(os.path.join(data_dir, 'Futuros.csv'),
parse_dates=['Fecha'],
index_col='Fecha',
date_parser=date_parser)
futuros = futuros.append(df.loc['2019-08-29'])
futuros.sort_index(inplace=True)
futuros['retorno'] = futuros['Cierre'].pct_change()
daily_returns_futuros = futuros['Cierre'].pct_change().dropna()
#Error en los retornos usando ARIMA
diff_arima = (daily_returns_futuros.reset_index()['Cierre'].values - mean_prediction[:13])
diff_abs_arima = np.abs(diff_arima)
np.sum(diff_abs_arima)/13
#Error en los retornos si se predice con el último valor
daily_returns_futuros.abs().sum()/13
#Error en el precio de cierre usando sólo ARIMA
predicted_close_arima = (1 + mean_prediction).cumprod() * df['Cierre'][-1]
diff_close_arima = (futuros.reset_index()['Cierre'][1:].values - predicted_close_arima[:13])
diff_close_arima_abs = np.abs(diff_close_arima)
np.sum(diff_close_arima_abs)/13
#Error en el precio de cierre si se predice con el último valor
diff_close_last_value = (futuros.reset_index()['Cierre'][1:].values - df['Cierre'][-1])
diff_close_last_value_abs = np.abs(diff_close_last_value)
np.sum(diff_close_last_value_abs)/13
Al aplicar ARIMA sobre la serie de retornos el resultado es prácticamente cero, lo cual es esperable ya que la idea es que el método capture la media condicional del proceso, y se ve claramente que ésta es (aproximadamente) cero en el gráfico.
En estas condiciones, la predicción a la que llegamos no es muy diferente a tomar el último valor de precio de cierre. De hecho, si comparamos el error absoluto medio (MAE) de ambas predicciones para las dos semanas posteriores a nuestros datos, la predicción que toma el último valor es un poco más precisa que ARIMA en este caso. Vemos también que GARCH nos da un valor de la volatilidad condicional estable, como buscábamos.
Últimos siete meses¶
Hacemos ahora las mismas predicciones pero restringiendo nuestros datos (la serie de tiempo de retornos) a los últimos siete meses. Esto se corresponde al período febrero-agosto, complementario a la siembra en Argentina (septiembre-enero).
df_seasonal_Feb_Aug = df_seasonal_Feb_Aug = df.loc[(df.index.month >1) & (df.index.month <9)]
df_seasonal_Feb_Aug['year'] = df_seasonal_Feb_Aug.index.year
g = sns.relplot(x='Fecha', y='Cierre', col='year', hue='year',
facet_kws=dict(sharey=False, sharex=False), col_wrap=3, legend=False,
kind='line', data=df_seasonal_Feb_Aug.reset_index())
for ax in g.axes:
plt.setp(ax.get_xticklabels(), rotation=45)
g.fig.suptitle("Cierre por año, ciclo Febrero-Agosto", size=16)
g.fig.subplots_adjust(top=.96, hspace=0.25);
retornos_seven_months = df.last('7M')['retorno'].dropna()
retornos_seven_months.plot();
model_arima_seven_months = sm.tsa.statespace.SARIMAX(retornos_seven_months.values,order=(p,0,q), seasonal=False)
fit_seven_months = model_arima_seven_months.fit()
mean_prediction_seven_months = fit_seven_months.forecast(n)
Nuevamente graficamos las predicciones de ARIMA de los retornos y el precio de cierre.
pd.Series(mean_prediction_seven_months).plot();
pd.Series((1 + mean_prediction_seven_months).cumprod() * df['Cierre'][-1]).plot();
residuos_seven_months = fit.resid
model_garch_seven_months = pf.GARCH(residuos_seven_months,r,s)
x = model_garch_seven_months.fit()
x.summary()
model_garch_seven_months.plot_predict(h=10)
model_garch_seven_months.plot_predict_is(h=50,figsize=(15,5))
#Error en los retornos usando ARIMA
diff_arima = (daily_returns_futuros.reset_index()['Cierre'].values - mean_prediction_seven_months[:13])
diff_abs_arima = np.abs(diff_arima)
np.sum(diff_abs_arima)/13
#Error en los retornos si se predice con el último valor
daily_returns_futuros.abs().sum()/13
#Error en el precio de cierre usando ARIMA
predicted_close_arima_seven_months = (1 + mean_prediction_seven_months).cumprod() * df['Cierre'][-1]
diff_close_arima = (futuros.reset_index()['Cierre'][1:].values - predicted_close_arima_seven_months[:13])
diff_close_arima_abs = np.abs(diff_close_arima)
np.sum(diff_close_arima_abs)/13
#Error en el precio de cierre si se predice con el último valor
diff_close_last_value = (futuros.reset_index()['Cierre'][1:].values - df['Cierre'][-1])
diff_close_last_value_abs = np.abs(diff_close_last_value)
np.sum(diff_close_last_value_abs)/13
Facebook Prophet¶
Entre los métodos alternativos que probamos, se encuentra la libreria Prophet de Facebook, que implementa un modelo de decomposición aditiva con tres componentes: trend ($g(t)$), estacionalidad ($s(t)$) y feriados (holidays, $h(t)$).
$$y(t) = g(t) + s(t) + h(t) + \epsilon_t$$Siendo $\epsilon_t$ un término de error que no ajusta el modelo. Este modelo es en esencia a un modelo aditivo generalizado que utiliza el tiempo como regresor.
Para más información, se puede consultar el paper original.
from fbprophet import Prophet
data = pd.DataFrame({'ds': df.index, 'y': df['Cierre']})
prophet_model = Prophet(weekly_seasonality=False)
prophet_model.add_seasonality(name='bianual', period=182, fourier_order=5)
prophet_model.fit(data);
prophet_forecast = prophet_model.predict(data)
prophet_model.plot_components(prophet_forecast);
ax = sns.lineplot(x=df.index, y=df['Cierre'], label='Precio de Cierre')
sns.lineplot(y=prophet_forecast['yhat'], x=prophet_forecast['ds'], label='Prophet forecast', ax=ax)
ax.set_xlabel('Fecha')
ax.set_ylabel('Precio')
ax.set_title('Cierre de Soja');
Calculamos a continuación, el AIC para el forecast de Prophet.
def AIC(ts, forecast, k):
sum_squared_errors = ((ts - forecast) ** 2).sum()
T = len(ts)
return T * np.log(sum_squared_errors / T) + 2 * (k + 2)
print('Prophet AIC =', AIC(df['Cierre'].values, prophet_forecast['yhat'].values, 3))
test_periods = pd.DataFrame({'ds': futuros.index})
prophet_predictions = prophet_model.predict(test_periods)
prophet_predictions.set_index('ds', inplace=True)
prophet_predictions['retorno'] = prophet_predictions['yhat'].pct_change()
last = df.iloc[-1]['Cierre']
ret = (prophet_predictions['yhat'][0] - last) / last
prophet_predictions['retorno'][0] = ret
prophet_mae = (prophet_predictions['retorno'] - futuros['retorno']).abs().sum() / (len(futuros) - 1)
print('Prophet forecast MAE =', prophet_mae)
ax = sns.lineplot(x=futuros.index, y=futuros['Cierre'], label='Precio de Cierre')
sns.lineplot(y=prophet_predictions['yhat'], x=prophet_predictions.index, label='Prediccion Prophet', ax=ax)
ax.set_xlabel('Fecha')
ax.set_ylabel('Precio')
ax.figure.autofmt_xdate()
ax.set_title('Cierre de Soja');
Autoregresiones Bayesianas¶
import pymc3 as pm
$AR(1)$ sobre precio de cierre¶
Probamos inicialmente modelando la serie de precios como un proceso $AR(1)$:
$$p_t = \phi_1 p_{t-1} + \epsilon_t$$sigma = 1.0
with pm.Model() as ar1:
phi1 = pm.Normal('phi_1', mu=0, sigma=sigma)
data = pm.AR('p', phi1, observed=df['Cierre'])
trace = pm.sample(10000, tune=4000)
map_ar1 = pm.find_MAP()
pm.traceplot(trace);
Como era esperable según lo visto en los plots de autocorrelación, el modelo asigna un valor cercano a $1$ para el coeficiente del primer término de lag.
A continuación, hacemos un plot de la predicción del modelo para los próximos 100 días.
phi1_hat = map_ar1['phi_1']
steps = 100
last_day = df.index[-1]
last_price = df.loc[last_day, 'Cierre']
forecast_period = pd.date_range(start=last_day + pd.DateOffset(1), periods=steps, freq='B')
forecast = np.repeat(phi1_hat, steps).cumprod()
forecasts = pd.DataFrame({'forecast': last_price * forecast}, index=forecast_period)
forecast_df = df.append(forecasts, sort=True)
def plot_forecast(forecast_df, steps=100):
ax = sns.lineplot(x=forecast_df.index[-2*steps:-steps], y=forecast_df['Cierre'][-2*steps:-steps], label='Precio de Cierre')
sns.lineplot(x=forecast_df.index[-steps:], y=forecast_df['forecast'][-steps:], label='Predicción', ax=ax)
ax.set_xlabel('Fecha')
ax.set_ylabel('Precio')
ax.figure.autofmt_xdate()
ax.set_title('Cierre de Soja');
plot_forecast(forecast_df)
$AR(2)$ sobre precio de cierre¶
Probamos ahora un modelo $AR(2)$, sabiendo que la autocorrelación parcial con los lags de mayor orden es cercana a 0.
$$p_t = \phi_1 p_{t-1} + \phi_2 p_{t-2} + \epsilon_t$$with pm.Model() as ar2:
phi1 = pm.Normal('phi_1', mu=0, sigma=sigma)
phi2 = pm.Normal('phi_2', mu=0, sigma=sigma)
data = pm.AR('p', [phi1, phi2], observed=df['Cierre'])
trace = pm.sample(10000, tune=4000)
map_ar2 = pm.find_MAP()
pm.traceplot(trace);
Vemos que el coeficiente del segundo término de lag es cercano a $0$, como era previsible. Veamos el plot de predicciones.
phi1_hat = map_ar2['phi_1']
phi2_hat = map_ar2['phi_2']
def ar2_process(steps, phis, last_prices):
prices = last_prices.copy()
for i in range(steps):
last, second_last, *rest = prices
next_price = phis[0] * last + phis[1] * second_last
prices = [next_price] + prices
return np.array(prices[:-2])[::-1]
second_to_last_price = df['Cierre'][-2]
forecast = ar2_process(steps, [phi1_hat, phi2_hat], [last_price, second_to_last_price])
forecasts = pd.DataFrame({'forecast': forecast},
index=forecast_period)
forecast_df = df.append(forecasts, sort=True)
plot_forecast(forecast_df)
$AR(1)$ sobre retornos¶
Probamos a continuación modelando los retornos como un proceso $AR(1)$
$$r_t = \phi_1 r_{t-1} + \epsilon_t$$with pm.Model() as ar1_ret:
phi1 = pm.Normal('phi_1', mu=0, sigma=sigma)
data = pm.AR('r', phi1, observed=df['retorno'].dropna())
trace = pm.sample(10000, tune=4000)
map_ar1_ret = pm.find_MAP()
pm.traceplot(trace);
Como era esperable, el coeficiente del primer lag de los retornos es aproximadamente $0$
phi1_hat_ret = map_ar1_ret['phi_1']
last_return = df['retorno'][-1]
forecast = (np.repeat(phi1_hat_ret, steps) * last_return + 1).cumprod()
forecasts = pd.DataFrame({'forecast': last_price * forecast}, index=forecast_period)
forecast_df = df.append(forecasts, sort=True)
plot_forecast(forecast_df)
Bayesian Structural Time Series¶
Otro enfoque más moderno en el análisis de series temporales es el que desarrolan Scott y Varian en el paper Predicting the Present with Bayesian Structural Time Series (2014).
A grandes rasgos, el modelo es un ensamble de predictores: combina la descomposición de series en diferentes variables de estado (tendencia, estacionalidad) con componentes regresivos.
Para la implementación, usamos Tensorflow Probability, una librería para programación probabilística creada sobre Tensorflow.
import tensorflow as tf
import tensorflow_probability as tfp
from tensorflow_probability import distributions as tfd
from tensorflow_probability import sts
tf.reset_default_graph()
trend = sts.LocalLinearTrend(observed_time_series=df['Cierre'])
seasonal = tfp.sts.Seasonal(num_seasons=33,
num_steps_per_season=123,
observed_time_series=df['Cierre'])
model = sts.Sum([trend, seasonal], observed_time_series=df['Cierre'])
with tf.variable_scope('sts_elbo', reuse=tf.AUTO_REUSE):
elbo_loss, variational_posteriors = tfp.sts.build_factored_variational_loss(
model, observed_time_series=df['Cierre'])
steps = 200
train_vi = tf.train.AdamOptimizer(0.1).minimize(elbo_loss)
with tf.Session() as sess:
sess.run(tf.global_variables_initializer())
for i in range(steps):
_, elbo_ = sess.run((train_vi, elbo_loss))
if i % 20 == 0:
print("step {} -ELBO {}".format(i, elbo_))
samples = sess.run({k: q.sample(50)
for k, q in variational_posteriors.items()})
for param in model.parameters:
print("{}: {} +- {}".format(param.name,
np.mean(samples[param.name], axis=0),
np.std(samples[param.name], axis=0)))
price_dist = tfp.sts.forecast(model,
observed_time_series=df['Cierre'],
parameter_samples=samples,
num_steps_forecast=steps)
num_samples = 10
with tf.Session() as sess:
price_mean, price_scale, price_forecast_samples = sess.run(
(price_dist.mean()[..., 0], price_dist.stddev()[..., 0],
price_dist.sample(num_samples)[..., 0]))
last_day = df.index[-1]
forecast_period = pd.date_range(start=last_day + pd.DateOffset(1), periods=steps, freq='B')
forecasts = pd.DataFrame({'forecast': price_mean}, index=forecast_period)
forecast_df = df.append(forecasts, sort=True)
colors = sns.color_palette()
c1, c2 = colors[0], colors[1]
ax = sns.lineplot(x=forecast_df.index, y=forecast_df['Cierre'], label='Precio de Cierre')
sns.lineplot(x=forecast_df.index, y=forecast_df['forecast'], ax=ax, label='Predicción BSTS')
ax.plot(forecast_period, price_forecast_samples.T, lw=1, color=c2, alpha=0.1);
ax.fill_between(forecast_period,
price_mean - 2 * price_scale,
price_mean + 2 * price_scale,
color=c2,
alpha=0.2)
ax.set_xlabel('Fecha')
ax.set_ylabel('Precio')
ax.figure.autofmt_xdate()
ax.set_title('Cierre de Soja');
Construimos a continuación un modelo que incorpora estacionalidad mensual y semanal a una tendencia lineal.
def build_model(ts):
estacionalidad_mensual = sts.Seasonal(num_seasons=12,
num_steps_per_season=246,
observed_time_series=ts,
name='estacionalidad_mensual')
estacionalidad_semanal = sts.Seasonal(num_seasons=52,
num_steps_per_season=7,
observed_time_series=ts,
name='estacionalidad_semanal')
tendencia = sts.LocalLinearTrend(observed_time_series=ts)
model = sts.Sum(
[estacionalidad_mensual, estacionalidad_semanal, tendencia],
observed_time_series=ts)
return model
tf.reset_default_graph()
model = build_model(df['Cierre'])
with tf.variable_scope('sts_elbo', reuse=tf.AUTO_REUSE):
elbo_loss, variational_posteriors = tfp.sts.build_factored_variational_loss(
model, df['Cierre'])
train_vi = tf.train.AdamOptimizer(0.1).minimize(elbo_loss)
steps = 200
with tf.Session() as sess:
sess.run(tf.global_variables_initializer())
for i in range(steps):
_, elbo_ = sess.run((train_vi, elbo_loss))
if i % 20 == 0:
print("step {} -ELBO {}".format(i, elbo_))
samples = sess.run(
{k: q.sample(50)
for k, q in variational_posteriors.items()})
price_dist = tfp.sts.forecast(model,
observed_time_series=df['Cierre'],
parameter_samples=samples,
num_steps_forecast=steps)
num_samples = 10
with tf.Session() as sess:
price_mean, price_scale, price_forecast_samples = sess.run(
(price_dist.mean()[..., 0], price_dist.stddev()[..., 0],
price_dist.sample(num_samples)[..., 0]))
Y ploteamos las predicciones con el nuevo modelo.
last_day = df.index[-1]
forecast_period = pd.date_range(start=last_day + pd.DateOffset(1), periods=steps, freq='B')
forecasts = pd.DataFrame({'forecast': price_mean}, index=forecast_period)
forecast_df = df.append(forecasts, sort=True)
colors = sns.color_palette()
c1, c2 = colors[0], colors[1]
ax = sns.lineplot(x=forecast_df.index, y=forecast_df['Cierre'], label='Precio de Cierre')
sns.lineplot(x=forecast_df.index, y=forecast_df['forecast'], ax=ax, label='Predicción BSTS')
ax.plot(forecast_period, price_forecast_samples.T, lw=1, color=c2, alpha=0.1);
ax.fill_between(forecast_period,
price_mean - 2 * price_scale,
price_mean + 2 * price_scale,
color=c2,
alpha=0.2)
ax.set_xlabel('Fecha')
ax.set_ylabel('Precio')
ax.figure.autofmt_xdate()
ax.set_title('Cierre de Soja');
La complejidad de este modelo hace que aquí lo estemos presentando como una caja negra. Aclarado esto, vemos que las predicciones de este modelo coinciden con los resultados vistos anteriormente: el mejor predictor del precio futuro es el último precio de cierre.
Modelo final¶
Predecir el comportamiento de los futuros de commodities ha sido motivo de estudio desde mediados del siglo pasado. En Commodity futures and market efficiency (2013, Kristoufek, Vosvrda), los autores señalan que los mercados de futuros de commodities (en particular, las relacionadas a la energía) muestran un alto grado de eficiencia, dada la multitud de actores que participan y hacen que el precio refleje toda la información disponible. Esto elimina cualquier tendencia observable en el precio.
Después de evaluar los distintos modelos observamos que todos coinciden en las predicciones: el mejor predictor del precio de cierrre en el período $t$ es el precio de cierre en el período $t-1$.
Para nuestra entrada a la competencia, nos inclinamos por un modelo ARIMA(1,0,1) sobre los retornos diarios. Siendo que este modelo es el mas simple e interpretable de todos los que pudimos evaluar, la navaja de Ockham nos llevo a elegirlo.
A continuación, construimos el dataset final para entrenar el modelo.
final_df = pd.concat([df, futuros.loc['2019-08-30':'2019-09-13']])
Entrenamos el ARIMA(1,0,1) sobre los retornos del dataset final.
arima_final = sm.tsa.statespace.SARIMAX(final_df['retorno'].values, order=(1,0,1), seasonal=False)
modelo_final = arima_final.fit()
Y finalmente creamos el csv con nuestras predicciones.
periodos = 20
forecast_retornos = modelo_final.forecast(periodos)
forecast_cierres = (1 + modelo_final.forecast(periodos)).cumprod() * final_df['Cierre'][-1]
forecast_fechas = pd.date_range(start=final_df.index[-1] + pd.DateOffset(1), periods=periodos, freq='B')
forecast_final = pd.DataFrame({
'Fecha': forecast_fechas.strftime('%d/%m/%Y'),
'Retorno': forecast_retornos,
'Cierre': forecast_cierres
})
forecast_final.to_csv(os.path.join(data_dir, 'forecast.csv'), header=False)
!cat data/forecast.csv